package com.platform.helper.maputil.geo;

/**
 * @所属包名称 :	com.diit.apppro.presentation.view.component.dialog
 * @描述 :   基于七参数坐标转换
 * @作者 :	kaki
 * @COPYRIGHT :	copyright(c) 2016,Rights Reserved
 * @版本 :	v1.0
 * @创建日期 :	2016/7/26
 */
public class CoordinateTransformSeven {
    /**
     * 转换过程中传递的三维坐标结果
     */
    private static double targetX = 0;
    private static double targetY = 0;
    private static double targetZ = 0;
    //最后转换返回三维坐标结果
    private static double lastX = 0;
    private static double lastY = 0;
    private static double lastZ = 0;
    /**
     * 高斯转换后传递的参数-经纬度
     */
    private double mTargetX = 0;
    private double mTargetY = 0;
    /**
     * 转换的参数-经纬度，高程
     */
    private double targetB = 0;
    private double targetL = 0;
    private double targetH = 0;
    //七参数
    private TransParaSeven transParaSeven = null;
    private boolean flg;
    /**
     * 80椭球和84椭球的长短半轴
     */
    final double Axis_XiAn80_a = 6378140.0;
    final double Axis_XiAn80_b = 6356755.28815753;

    final double Axis_WGS84_a = 6378137.0;
    final double Axis_WGS84_b = 6356752.3141999999;

    public CoordinateTransformSeven() {
        transParaSeven = new TransParaSeven();
        flg = false;
    }

    /**
     * 设置七参数、决定要不高思转换
     * @param para
     * @param isGauss
     */
    public void SetTransParaSeven(TransParaSeven para, Boolean isGauss) {
        transParaSeven = para;
        flg = isGauss;
    }

    /**
     * 基于七参数84坐标转换
     * @param longitude
     * @param latitude
     * @return
     */
    public CoordinateSeven TransformFromWgs84CoordinateBySeven(double longitude, double latitude) {
        CoordinateSeven coordinate = new CoordinateSeven();
        //经纬度的七参数转换，目标为西安80坐标系
        sevenParaTransformate(latitude, longitude, 0);
        //经纬度的高斯转换，得到平面直角坐标
        gaussBLtoXY(targetB, targetL);
        Coordinate.X = mTargetX;
        Coordinate.Y = mTargetY;
        return coordinate;
    }

    /**
     <summary>大地经纬度坐标转空间直角坐标</summary>
     <param name="sourceX">纬度，单位度</param>
     <param name="sourceY">经度，单位度</param>
     <param name="sourceZ">大地高</param>
     <param name="aAxis">参考椭球长半轴</param>
     <param name="bAxis">参考椭球短半轴</param>
     <param name="targetX">坐标X</param>
     <param name="targetY">坐标Y</param>
     <param name="targetZ">坐标Z</param>
    */
    /**
     * 第一步转换，大地坐标系换换成空间直角坐标系
     * @param B
     * @param L
     * @param H
     * @param aAxis
     * @param bAxis
     */
    private void BLHtoXYZ(double B, double L, double H, double aAxis, double bAxis) {
        double dblD2R = Math.PI / 180;
        double e1 = Math.sqrt(Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / aAxis;

        double N = aAxis / Math.sqrt(1.0 - Math.pow(e1, 2) * Math.pow(Math.sin(B * dblD2R), 2));
        targetX = (N + H) * Math.cos(B * dblD2R) * Math.cos(L * dblD2R);
        targetY = (N + H) * Math.cos(B * dblD2R) * Math.sin(L * dblD2R);
        targetZ = (N * (1.0 - Math.pow(e1, 2)) + H) * Math.sin(B * dblD2R);
    }

    /**
     * 第三步转换，空间直接坐标系转换为大地坐标系
     * @param X
     * @param Y
     * @param Z
     * @param aAxis
     * @param bAxis
     */
    private void XYZtoBLH(double X, double Y, double Z, double aAxis, double bAxis) {
        double e1 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(aAxis, 2);
        double e2 = (Math.pow(aAxis, 2) - Math.pow(bAxis, 2)) / Math.pow(bAxis, 2);

        double S = Math.sqrt(Math.pow(X, 2) + Math.pow(Y, 2));
        double cosL = X / S;
        double B = 0;
        double L = 0;

        L = Math.acos(cosL);
        L = Math.abs(L);

        double tanB = Z / S;
        B = Math.atan(tanB);
        double c = aAxis * aAxis / bAxis;
        double preB0 = 0.0;
        double ll = 0.0;
        double N = 0.0;
        //迭代计算纬度
        do {
            preB0 = B;
            ll = Math.pow(Math.cos(B), 2) * e2;
            N = c / Math.sqrt(1 + ll);

            tanB = (Z + N * e1 * Math.sin(B)) / S;
            B = Math.atan(tanB);
        }
        while (Math.abs(preB0 - B) >= 0.0000000001);

        ll = Math.pow(Math.cos(B), 2) * e2;
        N = c / Math.sqrt(1 + ll);

        targetZ = Z / Math.sin(B) - N * (1 - e1);
        targetB = B * 180 / Math.PI;
        targetL = L * 180 / Math.PI;
    }

    /**
     <summary>高斯转换时X是经度，Y是纬度</summary>
     <param name="x"></param>
     <param name="y"></param>
     <returns></returns>
    */
    /**
     * 第四部转换，高斯变换，大地坐标系转平面坐标系，84转80
     * @param mX
     * @param mY
     */
    private void gaussBLtoXY(double mX, double mY) {
        //参考椭球长半轴
        double m_aAxis = Axis_WGS84_a;
        //参考椭球短半轴
        double m_bAxis = Axis_WGS84_b;
        //中央子午线经度 济南117 威海123 巴州 87 带号
        double m_dbMidLongitude = TransParaSeven.daihao * 3;
        double m_xOffset = 500000;
        double m_yOffset = 0.0;
        try {
            //角度到弧度的系数
            double dblD2R = Math.PI / 180;
            //代表e的平方
            double e1 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_aAxis, 2);
            //代表e'的平方
            double e2 = (Math.pow(m_aAxis, 2) - Math.pow(m_bAxis, 2)) / Math.pow(m_bAxis, 2);
            //a0
            double a0 = m_aAxis * (1 - e1) * (1.0 + (3.0 / 4.0) * e1 + (45.0 / 64.0) * Math.pow(e1, 2) + (175.0 / 256.0) * Math.pow(e1, 3) + (11025.0 / 16384.0) * Math.pow(e1, 4));
            //a2
            double a2 = -0.5 * m_aAxis * (1 - e1) * (3.0 / 4 * e1 + 60.0 / 64 * Math.pow(e1, 2) + 525.0 / 512.0 * Math.pow(e1, 3) + 17640.0 / 16384.0 * Math.pow(e1, 4));
            //a4
            double a4 = 0.25 * m_aAxis * (1 - e1) * (15.0 / 64 * Math.pow(e1, 2) + 210.0 / 512.0 * Math.pow(e1, 3) + 8820.0 / 16384.0 * Math.pow(e1, 4));
            //a6
            double a6 = (-1.0 / 6.0) * m_aAxis * (1 - e1) * (35.0 / 512.0 * Math.pow(e1, 3) + 2520.0 / 16384.0 * Math.pow(e1, 4));
            //a8
            double a8 = 0.125 * m_aAxis * (1 - e1) * (315.0 / 16384.0 * Math.pow(e1, 4));
            //纬度转换为弧度表示
            //B
            double B = mX * dblD2R;
            //l
            double l = (mY - m_dbMidLongitude) * dblD2R;
            ////X
            double X = a0 * B + a2 * Math.sin(2.0 * B) + a4 * Math.sin(4.0 * B) + a6 * Math.sin(6.0 * B) + a8 * Math.sin(8.0 * B);
            //
            double ll = Math.pow(Math.cos(B), 2) * e2;
            double c = m_aAxis * m_aAxis / m_bAxis;
            //N
            double N = c / Math.sqrt(1 + ll);
            //t
            double t = Math.tan(B);
            double p = Math.cos(B) * l;
            double dby = X + N * t * (1 + ((5.0 - t * t + (9.0 + 4.0 * ll) * ll) + ((61.0 + (t * t - 58.0) * t * t + (9.0 - 11.0 * t * t) * 30.0 * ll) + (1385.0 + (-31111.0 + (543 - t * t) * t * t) * t * t) * p * p / 56.0) * p * p / 30.0) * p * p / 12.0) * p * p / 2.0;
            double dbx;
            dbx = N * (1.0 + ((1.0 - t * t + ll) + ((5.0 + t * t * (t * t - 18.0 - 58.0 * ll) + 14 * ll) + (61.0 + (-479.0 + (179.0 - t * t) * t * t) * t * t) * p * p / 42.0) * p * p / 20.0) * p * p / 6.0) * p;
            mTargetX = dbx + m_xOffset;
            mTargetY = dby + m_yOffset;
        } catch (Exception ex) {
        }
    }

    /**
     * 七参数转换，包含了第二步转换，空间直角坐标系里七参数转换
     * @param B
     * @param L
     * @param H
     */
    private void sevenParaTransformate(double B, double L, double H) {
        double m_sourceAxis_A = Axis_WGS84_a;
        double m_sourceAxis_B = Axis_WGS84_b;
        double m_targetAxis_A = Axis_XiAn80_a;
        double m_targetAxis_B = Axis_XiAn80_b;
        try {
            //根据坐标类型获取需要的空间直角坐标系XYZ
            //double X = 0.0, Y = 0.0, Z = 0.0;
            BLHtoXYZ(B, L, H, m_sourceAxis_A, m_sourceAxis_B);
            //7参数计算
            double Ex, Ey, Ez;


            Ex = TransParaSeven.m_dbXRotate / 180 * Math.PI;
            Ey = TransParaSeven.m_dbYRotate / 180 * Math.PI;
            Ez = TransParaSeven.m_dbZRotate / 180 * Math.PI;

            double newX = TransParaSeven.m_dbXOffset + TransParaSeven.m_dbScale * targetX + targetY * Ez - targetZ * Ey + targetX;
            double newY = TransParaSeven.m_dbYOffset + TransParaSeven.m_dbScale * targetY - targetX * Ez + targetZ * Ex + targetY;
            double newZ = TransParaSeven.m_dbZOffset + TransParaSeven.m_dbScale * targetZ + targetX * Ey - targetY * Ex + targetZ;

            //根据坐标类型获取需要的空间直角坐标系XYZ
            lastX = newX;
            lastY = newY;
            lastZ = newZ;
            XYZtoBLH(newX, newY, newZ, m_targetAxis_A, m_targetAxis_B);
        } catch (Exception ex) {

        }
    }

}
